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ABSTRACT 

A scheme is presented for accurately propagating the gravitational field constraints 
in finite difference implementations of numerical relativity. The method is based on 
similar techniques used in astrophysical magnetohydrodynamics and engineering elec- 
tromagnetics, and has properties of a finite differential calculus on a four-dimensional 
manifold. It is motivated by the arguments that 1) an evolutionary scheme that nat- 
urally satisfies the Bianchi identities will propagate the constraints, and 2) methods 
in which temporal and spatial derivatives commute will satisfy the Bianchi identities 
implicitly. The proposed algorithm exactly propagates the constraints in a local Rie- 
mann normal coordinate system; i.e., all terms in the Bianchi identities (which all vary 
as d^g) cancel to machine roundoff accuracy at each time step. In a general coordi- 
nate basis, these terms, and those that vary as dg d^g, also can be made to cancel, but 
differences of connection terms, proportional to (dg)^, will remain, resulting in a net 
truncation error. Detailed and complex numerical experiments with four-dimensional 
staggered grids will be needed to completely examine the stability and convergence 
properties of this method. 

If such techniques are successful for finite difference implementations of numer- 
ical relativity, other implementations, such as finite element (and eventually pseudo- 
spectral) techniques, might benefit from schemes that use four-dimensional grids and 
that have temporal and spatial derivatives that commute. 

Subject headings: relativity: numerical — black holes 



1. Introduction 

The quest for solutions of dynamical strong gravity problems, such as black hole formation 
or mergers and gamma-ray burst production, that are astrophysically relevant and accurate enough 
to predict gravitational wave forms, has occupied much of the last half of the 20* century and 
beyond. Its history has been frustrated by the need to address several unforeseen numerical prob- 
lems, including 1) proper initial data to begin the evolution, 2) the development of coordinate 
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and physical singularities during the evolution, and 3) the growth of numerical instabilities in the 
time-dependent solutions. These problems have been dealt with by addressing each in turn with 
such techniques as 1) elevation of the initial data problem to a complete sub-field; 2) develop- 
ment and analysis of appropriate gauge/slicing conditions that avoid coordinate singularities; 3) 
excision of black hole centers, inside horizons, to avoid physical singularities; and 4) the use of 
symmetric/hyperbolic equations to enhance numerical stability. 

One of the several remaining problems in this field is that, for some problems and some co- 
ordinate systems, the constraint-violating modes will grow exponentially. These eventually over- 
whelm any solution in only a short period of time (< 100 M), rendering a long simulation (and 
the computation of any gravitational wave forms) impossible. For very high-resolution simula- 
tions the exponential growth of errors begins early and continues until the errors diverge. On the 
other hand, for simulations with coarse spatial resolution the errors begin and remain at the trun- 
cation level until the much smaller constraint violations grow to a level that exceeds the truncation 
accuracy. Then the solution joins the general exponential growth seen in the high-resolution sim- 
ulations, blowing up in much the same manner as in the high-resolution case (Scheel et al. 2002). 
The similar behavior of this error at a variety of mesh spacings indicates that it may be a numerical 
solution to the discrete equations that are being integrated. Current attempts to solve this problem 
include adding the constraints as penalty functions to the evolution equations and techniques that 
re-converge the constraint equations every few time steps. 

Constraint propagation is also an issue in the solution of Maxwell's equations. Techniques 
for doing so in the fields of astrophysical magnetohydrodynamics (MHD) and electromagnetics of 
antennas and waves have been in place for decades. These enforce the constraints not just stably, 
but to machine accuracy. Finite difference methods for constraint propagation in astrophysical 
MHD are known as the Evans-Hawley constrained transport (CT) method (Evans & Hawley 1988) 
and now are an integral part of publicly used codes, such as ZEUS (Stone & Norman 1992a; 
Stone & Norman 1992b) and ZEUS-3D (Clarke 1996). In engineering electromagnetics these are 
known as the Yee algorithm (Yee 1966), and have many variants (De Raedt et al. 2002). They all 
involve building a mesh that is staggered in space, and often in time as well, and then defining 
appropriate vector and scalar quantities at whole or half mesh points. 

While the success of CT for electromagnetics is certainly due in part to the linearity of the 
physical equations, its ability to maintain accuracy of the solenoidal (V • 5 = 0) and Coulomb 
(V ■£ = 4;tpt) constraints, to a few parts in 10^^^^^ over tens or hundreds of thousands of time 
steps, is enticing. If such an algorithm can be found for numerical relativity, it could be as useful 
as excision and other such proven methods. 

In this paper the properties of the electromagnetic CT method (and of spacetime itself) are 
examined, and a similar method is developed for numerical relativity. It is the thesis of this paper. 
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and subsequent ones in this series, that CT methods work because they build spacetimes in which 
the temporal and spatial derivatives commute ([3o,3!] = 0(8^), where £^ is the machine roundoff 
error). This naturally enforces the Bianchi identities, and it is those implicit identities that prop- 
agate the constraints. This thesis is not fully tested in this first paper on the subject. In fact, it 
will take some time and numerical effort to verify or refute it. Instead, only the first steps are 
taken here. It is shown that the proposed CT scheme for NR works exactly in Riemann normal 
coordinates; in general coordinates most terms also cancel. Detailed numerical simulations in four 
dimensions will be needed to fully explore the method's stability properties to see if these condi- 
tions are sufficient for stable constraint propagation. If successful, however, similar CT methods 
should be possible for other implementations of numerical relativity, not just for finite differences. 

This paper is intended to be the first in a long series that will culminate in a numerical code 
that is capable of simulating black hole formation and the gamma-ray burst jet generation and 
gravitational wave production that is expected to result from such events. These issues are im- 
portant both for dealing properly with the energetics in the system as well as with the expected 
observational consequences of the event. In order to treat these problems properly, such a code 
must be capable of evolving the relativistic gravitational field, as well as the fluid matter flowing 
within that field and the electromagnetic field generated by currents flowing in that matter. Part of 
achieving this goal will be to present a consistent numerical method, from the relativistic gravita- 
tional field to specific issues of stellar mergers and collapse. In this paper we lay the groundwork 
for generating the time-dependent gravitational field and for evolving the electromagnetic field in 
that metric. In subsequent papers we will present tests of the constrained transport techniques de- 
veloped herein, and eventually add the stress-energy due to matter and fluid motion to complete 
the code. Then, specific astrophysical problems will be addressed. The ultimate aim of this work 
is to foster a closer relationship between astronomers who observe black hole systems and those 
numerical physicists and astrophysicists who study them theoretically. 

2. Review of CT for Electromagnetics in Flat Spacetime 
2.1. The Evans-Hawley CT Method for MHD 

For astrophysical MHD the field equations that are solved are 

B = -cVxE (1) 

with the solenoidal constraint 

V-fi = (2) 
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(Additional equations are solved, of course, including the conservation of mass, momentum, and 
energy; but these are not relevant here.) The Evans-Hawley constrained transport technique sat- 
isfies equation (2) on the initial hypersurface, usually to machine accuracy, and then uses a dif- 
ferencing scheme for equation (1) that ensures that equation (2) is satisfied on each subsequent 
hypersurface to the same level of accuracy as on the first. This is done by staggering the grid in 
space and time (see Figure 1). At whole time steps the magnetic field vector components are de- 
fined normal to, and centered on, grid cube faces. At half time steps the electric vector is defined 
parallel to, and centered on, cube edges. For the initial conditions the magnetic field is derived 
from a vector potential 

B = VxA (3) 

This vector potential is defined on cube edges at t = ^t, and for evolutionary computations B = 
dB/dt is defined on cube faces, centered in time between "B and i.e., at t = "'^^t. 

At t— we see that the solenoidal constraint is satisfied to machine accuracy by this method. 
For Ax = Ay = Az we have, simply, 

VB = +Bx — -Bx + +By — -By + +B^ — -B^ 
= V-VxA 

= (++^z - +-^z) - i++Ay - +-Ay) - i-+Az - —Az) + - — ^y) + 

i++Ax - - i++Az - -+Az) - i-+Ax - —Ax) + i+-A^ - — + 

{++Ay - -+Ay) - {++Ax - -+Ax) - i+-Ay - —Ay) + {+-Ax - —Ax) 

= 0(8,) (4) 

where +5, and Hi are components on upper and lower cube /-faces, respectively, and the two pre- 
appended subscript signs on vector potential components Ay indicate the 7* edge at the intersection 
of the upper and/or lower cube faces with normals in the two spatial dimensions orthogonal to j. 
Similarly, taking the numerical divergence of equation (1) at 2?, we have 

V B = -cV-VxE = 0(8,) 

so, at time t = ^tthe magnetic field remains divergence-free 

V^B = V- °B + A? 2(V-fi) =0(8,) 



because it is the sum of two divergence-free fields. The solenoidal constraint is therefore preserved 
to machine accuracy specifically because the vector identity V -V x E = is naturally satisfied by 
the differencing scheme. 
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In the above equations the electric field can be computed by any means, and CT still would be 
maintained. Ideal MHD codes use Ohm's law with infinite conductivity to relate the fluid velocity 
and the magnetic field: 

V 

E = — xfi 

c 

Some astrophysical codes, like ZEUS and several in Japan, use the method of characteristics to 
preserve Alfven waves in the E update, and so are often called MOC-CT codes. However, it is 
important to realize that CT results simply from the grid staggering and has nothing to do with the 
method of characteristics or any other E update scheme. 



2.2. The Yee Algorithm for Electromagnetics 

In full electrodynamics two more Maxwell's equations are used to determine the electric field, 
instead of Ohm's law. In vacuum these are 

E = cVxB - AtU (5) 
VE = 4%pc (6) 

Figure 1 then must be modified to include the new quantities J (current density), Pc (charge den- 
sity), and (|) (electric potential). (See Figure 2.) We also add another initial data problem at t = ^2t 
that solves equation (6) for E. If, for the moment, we choose the Coulomb gauge and ignore the 
vector potential, we have 

V^^ = -47rpc 
E = -V(^ 

Electric potential and charge density then must be defined at cube comers on half time steps, and 
current density must be defined on cube edges at whole time steps (same as curlB, A, and E; see 
Figure 2). In an arbitrary gauge, we will have 

E = -V(^ - -A (7) 
c 

indicating that A and E are co-located. The equations are closed by specifying the evolution of Pc 
and J, which together must satisfy the conservation of charge 

V-7 = -Pc (8) 

so pc must be defined on cell comers at whole time steps. This method of staggering was suggested 
buy Yee almost 40 years ago (Yee 1966). 
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Fig. 1. — Space-time representation for the Evans-Hawley CT scheme. Open circles are face- 
centered on the cubes and crosses are edge-centered. 




Fig. 2. — Space-time representation for the Yee algorithm. Similar to Figure 1, but with filled 
circles located at cube corners. Another time step has been added C^t, along with electric field 
quantities. 
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Constraints are preserved in the Yee algorithm in the same manner as in the Evans-Hawley 
algorithm. For Faraday's law and the solenoidal constraint, the procedure is identical. And for 
Ampere's law we have 

VE = c V- VxB - 47rV-7 

= 0(£,) + 

which gives the following update for divE: 

3(V-£) = -^{V-E) + \V-E) At 

= An ^p,. + 0{Zr) 

So the Coulomb constraint is preserved to machine accuracy as long as is conserved in the 
update. 



2.3. Covariant formulation of CT for Electrodynamics 

Figure 3 re-casts the Yee algorithm in covariant form, using the Faraday tensor F (instead 
of the vector fields), the four-current J, and the vector four-potential A. This will give important 
clues to developing a CT scheme for Einstein's field equations. Maxwell's equations then become, 
including constraints. 



V-M 
VF 





47rJ 



(9) 
(10) 



where M = *F is the Maxwell tensor (the dual of F) and V is now the /oMr-gradient operator 

3 3 3 9 



V = (^,J^,^,^). Because M and F are antisymmetric, they satisfy tensor identities (analogous 
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Foi ' ji, /ti ' Jo Foi Jo 

Fig. 3. — Same as Figure 2, but with covariant notation. Note placement of vector components (on 
hypercube edges) and Faraday tensor components (on hypercube faces). 
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to the vector identities V -V xE) called the Bianchi identities 

V-(V-M) = 
V-(V-F) = 

These identities are related to J. A. Wheeler's classic statement that "the boundary of a boundary 
is zero". But, does the staggered grid in Figure 3 automatically satisfy these identities to machine 
accuracy? A quick analysis of V ■ F shows that, in fact, it does. V ■ F is a vector that is defined 
on cell edges at whole time steps and on cell corners at half time steps, and it involves tensor 
components that are two half-steps away from scalar points (whole time-step cell comers). Taking 
the divergence of this vector causes like components to cancel, so that 

V-(V-F) = 0{er) 

holds in this differencing scheme. 

Because of the zero on the right hand side of equation (9), that equation itself also is often 
described as a Bianchi identity 

dF = (11) 
where d¥ is the differential of the tensor F, which in component form is given by 

(^)aPY = P[ap,Y] = Fap,Y + Fpy," + ^ya,p 

where the comma denotes ordinary differentiation (Fc^p y = dFf^^/dx^) and the brackets denote 
permuted summation.^ This allows the Faraday tensor to be derived from a vector four-potential 

¥ = dA 

or 

This is the covariant form for equations (3) and (7). Does the staggered grid automatically enforce 
ddA — also? Yes. We have already shown this to be the case for the magnetic part (equation 4). 
For the electric part of ddA = we have 

VxE = -V X V(|) - V X -A 

c 

= -B + 0{er) (12) 
c 



'In this paper Greek indices range from to 3, Roman indices i, j, k, ... range from 1 to 3, and Roman indices 
a, b, c, ... will be used to denote a set of three integers with one of the spatial indices missing {i.e., one of the sets [0, 
1, 2], [0, 1, 3], or [0, 2, 3]). 
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the second term is B/c to order and the first term is zero to the same order because of the 
staggered grid (see the t=~^t panel in Figure 2). So equation (12) is just the same Faraday's law 
that we are solving to machine accuracy in the spatial part of equation (9) (or equation 1). 

To summarize, then, the staggered grid naturally satisfies the Bianchi (and vector) identities 
in space and time 

F"^pa = (13) 
(Aoc,p - Ap^„),y + (Ay,cx - Acx,y),p + (Ap,y - Ay^p),a = (14) 

to machine accuracy for any antisymmetric tensor F and for any four-vector A. Here we use 
the Einstein summation convention, where a repeated index indicates summation over the four 
coordinates (F"P p = Ep^Q3F°'P/3xP). Note that equation (13) uses the raised form of F, but in 
flat space this involves only multiplying by ±1 with the Minkowski metric. As a result of this 
cancellation, when the spatial parts of equations (9) and (10) are integrated forward in time 

p.(V-F) = 47IP-J P-(V-M) = 

(where the spatial projection tensor P'^P — n'^nP + g'^P is orthogonal to n), the time parts (the 
constraints) are automatically satisfied to machine accuracy 

n-(V-F) = 47W1-J n-(V-M) = 

with no additional computation required. A staggered grid, therefore, has "deep geometric signifi- 
cance", because it naturally satisfies the Bianchi identities. 



3. CT for Electrodynamics in Curvilinear and Curved Spacedme 

CT also works in curvilinear coordinates and in curved spacetime, but the method requires 
iteration. Again, Maxwell's equations are 

V-F = 47tJ d¥ = Q 

but the gradient operator is now the covariant derivative rather than the ordinary derivative 

F«P.p = 47iJ« F[„p.y] = (15) 

where 
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are the covariant derivatives of a general second-rank tensor T, 

r\ = g«^r^p^ (16) 

TaPY = 2 (SctP'Y + SPa.Y^gPY.a) ^^^^ 

are different forms of the connection coefficients, and g"P is the inverse of the metric g^p. However, 
because the Faraday tensor is antisymmetric (as is the Maxwell tensor), equations (15) reduce to 

F'"% = 47iJ'« F[„p,,] = (18) 

where 




and ^/—g is the volume element (square root of the negative metric determinant). Equations (18) 
involve only simple differences and known values for the metric. However, unless g is a diag- 
onal, the raised version of the Faraday tensor F'^'P involves the sum of several lowered version 
components at different grid locations 

pap = g«^gPvF^^ 

Therefore, the time update of the fundamental variables F^v will necessarily be implicit, and there- 
fore iterative, as it involves sums over time as well as space. 

It is important to realize, however, that even in curvilinear coordinates and curved spacetime, 
the Bianchi identities 

F'"%a = (19) 

will be satisfied to machine accuracy, because F' is constructed before it is differenced and 
because that differencing is done in precisely the same manner as in equation (13). It does not 
matter that F' involves the sum of many tensor components and products of metric components. 
It only matters that the |3 and a derivatives commute. 



4. A General Finite Difference Prescription for CT for Tensor Field Evolution Problems 

Figure 3 suggests the following geometric prescription for a staggered grid when solving 
covariant tensor field evolution problems. This prescription is depicted schematically in Figure 4, 
and appears to work for tensors up to at least rank five. The basic rules are 
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n+U 



Scalars: 

(hyper cube 
corners) 



Vectors: 

(hypercube 
edges) 



2-Tensors: 

(hypercube faces 
& corners) 



3-Tensors, 
Vg, r, etc. : 

(hypercube edges & 
cube body centers) 



4-Tensors: 

(hypercube corners, 
faces, 

& body centers) 




Fig. 4. — Graphical depiction of the four-dimensional staggered grid scheme for different tensor 
quantities. (See Section 4 for a full description.) The first panel shows filled circles where scalars 
are located (hypercube comers.) Heavy arrows in the second panel show how to locate the fol- 
lowing vector quantities: Jo (shift one-half step in the dimension); Ji (shift one-half step in the 
1 dimension). The third panel shows how to locate R03 by shifting one-half step in the and 3 
dimensions, and R13 by shifting in the 1 and 3 dimensions. Remaining panels demonstrate third 
and fourth ranked tensors. Note the hypercube body-center location of Rq 1 23 • 
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1. Extend and stagger the grid in time as well as in space. The time extension need not be very 
deep — only enough cells to compute the tensor components, derivatives, etc. to the order 
of the method. In this paper we use a second order differencing scheme, so we need only 
one additional half, and one full, time slice. 

2. Treat time otherwise like a spatial coordinate. That is, use the same differencing scheme in 
time as used in space so that temporal and spatial differences commute. 

3. Starting at comer nodes on the four-dimensional hypercube cell, define the following quanti- 
ties in the manner described below. Simply put, a tensor component is located one-half step 
away from hypercube cell comers in each dimension specified by that component's indices; 
two repeated indices are equivalent to no shift at all. 

(a) Scalars: located on 4-cube cell comers (3-cube corners at whole time steps) 

(b) Vectors: located on 4-cube edges, shifted one-half cell step in the direction specified 
by that component. That is, 

i. j'^ on 3-cube cell comers at half time steps 

ii. centered on 3-cube x edges at whole time steps 

iii. centered on y edges 

iv. centered on z edges 

(c) One-forms: defined like vectors 

(d) Second-ranked tensors of any type 

i. Raa: defined like scalars, on 3-cube cell comers at whole time steps 

ii. Ro;: centered on 3-cube /-edges at /za//time steps 

iii. R(y: centered on ij faces at whole time steps 

(e) Third-ranked tensors and connection coefficients 

i- Taap, T^^^, T^a.a- defined like JP 

ii. Foij, T/oj, r,jo: centered on 3-cube «j-faces at half time steps 

iii. Tijk'. 3 -cube-centered at whole time steps 

(f) Fourth-ranked tensors 

i. Raapy- defined like Rpy 

ii. RapY5 (« 7^ P 7^ Yt^ 5): located at 4-cube body centers 

(g) Bianchi identities: defined like third-ranked tensors (at least one index must be re- 
peated). Examples are Roo23;2 (same as J^) and Roi23;2, (centered on the 1-3 face at 
half time steps). 
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This differencing scheme has a number of properties that make it look like a finite implemen- 
tation of differential calculus. First, the differential operator, which creates a tensor of one higher 
order (e.g., equation 11), naturally places the new tensor on the proper grid if the differencing is 
centered. This is also true of the generation of T^^^y from the metric field g^p. Second, one of 
the most fundamental properties of a spacetime, the covariant derivative of the metric (gap;Y) 
ishes, is naturally satisfied to machine accuracy because of this property. Third, the contraction of 
a mixed tensor is trivial: for each staggered component of the contracted tensor the four compo- 
nents of the parent tensor that are needed for the sum are already located at the same grid point as 
the contracted component. No additional averaging is needed. Fourth, as shown below, in a local 
Riemann normal coordinate system, the Bianchi identities are satisfied to machine accuracy. All 
terms cancel exactly, so the constraints are propagated exactly as well. 

Certain special tensors also have interesting properties. The Kronecker delta 5"p, for exam- 
ple, has non-zero elements (unity) only at hypercube cell comers. This is also true of other identity 
tensors (S^^^/f S"P\^V' etc.), which are ±1 at cell comers. The Levi-Civita tensor e^^y^ ^nd anti- 
symmetric symbol [ocPyS] are just the opposite. They are zero everywhere except at the hypercube 
body centers. They look much like the identity tensors, but on a grid that is shifted one-half step in 
each dimension. Furthermore, as the Levi-Civita tensor expresses the volume element 

its placement at the hypercube center creates a natural scheme for forming volume integrals over 
those hypercubes. 



5. CT for Numerical Relativity 

5.1. Statement of the Problem 

For reasons that are developed more fully below, we will use mixed tensors to define the 
problem of numerical relativity. As noted above, such tensors will lend themselves easily to con- 
traction. 

The classic problem of general relativity is to solve Einstein's field equations 

G"p = 87tT«p (20) 

(c = G = 1) for the metric coefficients. The Einstein curvature tensor is derived from the Ricci 
curvature tensor 



G«p ^R«p-i5«pi? 



(21) 
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with 

R = R^'a (22) 

being the Ricci curvature scalar. The Ricci tensor is the contraction of the Riemann tensor on the 
first and third indices 

R«p = R^%p (23) 

The Riemann tensor is the full statement of curvature of the spacetime. In its mixed form it is given 
by 

R% = r«\,-r«^6 + r%rP^-r%rP,5 (24) 

for a coordinate basis. The doubly-raised connection coefficients are given by 

pap^ = gP^r% (25) 

We assume that the metric g^^ has a unique inverse g"P such that 

g^'^g^p = 5«p (26) 

In this paper we will treat only the vacuum problem (the source of stress-energy T°'p = 0) so that 
equation (20) becomes 

R«p = G"p = (27) 

The Riemann tensor R^l^ys possesses several symmetries, including algebraic antisymmetry on a 
and P (and on y and 5) and differential symmetries (Bianchi identities) 

R%8,,^ ^ R«P,5;e + R"Pey;5 + R"p8e;y = (28) 

The reader will note that the mixed Riemann tensor is missing one additional symmetry that is 
possessed by the covariant version: Rapyg = RySap- (A raised index cannot be swapped with a 
lower one.) When contracted on the first and third indices, the Bianchi identities become simply 

R"p;a - ^^,p = G«p.„ = (29) 

i.e., the divergence free condition on the Einstein tensor. These four conditions are responsible for 
propagating the four constraints 

G^p = (30) 

if the other equations are satisfied 

G'p = (31) 
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5.2. CT in Riemann Normal Coordinates 

At any point P in spacetime one can construct many transformations L"^ to locally Lorentz 
systems such that 

gdp = L^'&L^pg^v = Tl^p (32) 

in a neighborhood of that point. However, only one of those systems — the Riemann normal sys- 
tem — also has vanishing gradients of the metric and, therefore, vanishing connection coefficients 
in that same neighborhood 

= 

In this coordinate system, in the neighborhood of P, covariant derivatives become ordinary deriva- 
tives and the Riemann tensor and its Bianchi identities become 

R^s = r^^5,t-r\5 (33) 
R^p.„. „ = r«Ps - r'^P. + T^K - T^K + r«P. - t^^, = o (34) 

[Yo;e] 5, ye y,oz y,eo e,yo £,07 o,ey ^ ^ 

In the proposed staggered grid scheme in the previous section, each of the terms of R^'^j^g.g] would 
be evaluated at the same grid point, because they each have the same five indices. So the sum 
can be accomplished without additional averaging from other grid points. We further note that 
each term has a duplicate with the opposite sign, differing only in the order of the derivatives {e.g., 
^"'^5 ye ~ ^"'^5 ey-*' '^^^ therefore draw the following conclusion: If a numerical scheme is 
constructed such that derivatives commute (both space-space and space-time), then in Riemann 
normal coordinates the equation for the Bianchi identities (34) will be satisfied to machine accu- 
racy, resulting in the propagation of the constraints to machine accuracy. We note that the scheme 
proposed in Section 4 possesses the required properties. 

How does satisfying the Bianchi identities propagate the constraints numerically? This is easy 
to show in Riemann normal coordinates. For the vacuum problem, the constraints are given by 

ROp = (35) 

and we seek a scheme in which the constraints propagate to machine accuracy 

r\5 = 6>(£,) (36) 

But satisfying the Bianchi identities (equation 34) will mean that the contracted identities are also 
satisfied to machine accuracy 

R^p,a = 0{Er) (37) 
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or 

RVo = -^^J+Oier) (38) 

All that needs to be shown is that R'a ^ = 0(er) also. For the three momentum constraints 

(P = 7) this is straightforward, because R'y = are the field equations being computed, so R'j 

is 0{er) by definition.^ So the spatial gradients R'j f will be O(e^), thereby propagating R^j to 
machine accuracy. 

For the Hamiltonian constraint (P = 0) we require satisfaction of the momentum constraints 
on the hypersurface, i.e., R'^ •> = 0{Zr)- Therefore, as long as the momentum constraints propagate 
to machine accuracy on each hypersurface, which we have shown above to be the case, the Hamil- 
tonian constraint will propagate also. Note, however, if the gradient of the momentum constraints 
has a constant bias, then the Hamiltonian constraint will grow. On the other hand, if the divergence 
fluctuates randomly, the the Hamiltonian constraint will fluctuate only randomly as well. 

5.3. CT in a General Coordinate Basis 

Because a Riemann normal coordinate system is local only and cannot be used to cover the 
entire spacetime, we are forced to deal with non-zero connection coefficients. But we still will 
attempt to propagate the constraints in the same manner — with a scheme that enforces the Bianchi 
identities by employing temporal and spatial differences that commute. 

5.3.1. Propagation of the Constraints 

It is important to note that the Bianchi identities are never explicitly calculated in the CT 
method. Instead, we develop a scheme in which 

^This statement requires a little clarification. Of course, the solution of R'j = is accurate to only 0{Ztr)- 
However, if we use this truncation-accurate solution to re-compute R' j, using exactly the same mathematical definition 
of terms that we used in the evolution equation, then that re-computed R' j will be zero to machine accuracy. (It will not 
be so only if we use a different differencing scheme than the one used in the original evolution equation.) As a simple 
example, consider a hue of code that computes y = ax + b. Then, if we later compute the function f = y — ax — b, 
by definition, / will be zero to machine accuracy. 
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is satisfied to at least truncation accuracy for all time. (We will find that machine accuracy may 
be an unattainable goal.) However, in a general coordinate system, the update of R°p will depend 
on an "advection" of curvature from surrounding cells, due to the covariant derivative connection 
terms. So, we equivalently seek a scheme in which 

ROp.o = (39) 

to at least truncation accuracy. If we set up a grid in which 

R"^Y6;e] = 0(£,,) (40) 

then we can proceed in much the same manner as in Riemann normal coordinates, but with covari- 
ant derivatives. If the Bianchi identities are satisfied to truncation accuracy, then their contracted 
form also will hold 

R%;0 = -R^;; + 0(e,,) 
R^-;0 = -^U-j + 0{etr) 

So, if the momentum constraints are satisfied on each hypersurface, the Hamiltonian constraint 
will propagate, albeit along geodesies, not coordinate lines. And the momentum constraints will 
propagate if the field equations are satisfied to at least truncation accuracy. 

5.3.2. Cancellation of the d^T Terms 

In order to demonstrate cancellation of terms in the Bianchi identities, we choose the follow- 
ing form for the mixed Riemann tensor 

R"^5 = ^ (r«p5,y - r«Py,5 - rP\, + rP\5 

+ r^/^grP^T - r%rP^5 + ^.^r^'j - r«^rP^) (4i) 

with the singly-raised connection coefficient re-computed from the doubly-raised ones and the 
gradient of the inverse metric 

rP,, = ^(rP\-rP,-gPv,,) (42) 

This form has the following properties 

1. R*^Py5 possesses explicitly all of the algebraic symmetries discussed in Section 5.1. 
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2. When the Bianchi identities are formed^ 

R^^Sej + r^yR^Pge + rP^yR^e = (43) 

and equation (42) is inserted, the F^'^y^se terms will cancel to machine accuracy, as described 
before in section 5.2. 



5.3.3. Cancellation of the TdF Terms 

In equation (43) the r'^^eP^'^g y terms will cancel explicitly algebraically. Will they cancel 
numerically also? The answer is yes, if we apply the following numerical procedures 



1. Use linear interpolation (averaging) to determine quantities at intermediate grid points 

2. When forming a product, such as g^v T^^^y, first average the factors to the grid point in ques- 
tion, then form the products and finally the sum. Do not form the products on different grid 
points and then average. 

Consider, for example, the two following Bianchi identity terms 

r"^5,erP^y - r^yr«^5,e (44) 

Algebraically, of course, the two terms cancel. However, in a staggered grid scheme they are not 
computed in the same manner. The first term comes from differencing a FT term in the first term 
of equation (43), while the second comes from the connection of one of the dT terms in the last 
term of that same equation. One is the difference of an average, while the other is the average 
of differences. But, with linear averaging, we see that, under these conditions, the chain rule is 
satisfied to machine accuracy 




^Note that, because of the antisymmetry in the last two indices of R^l^yg and the symmetry in the last two indices of 
r"p y, the connection terms involving y and 5 wiU cancel, just as they did for the antisymmetric Faraday tensor Bianchi 
identities in equation (18). 
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+ 



(r«^,erP,,) + (r«'^5rP,,,e) 



(45) 



(46) 



where the pre-appended subscript ^ signifies the spatial node position at which that quantity is 
evaluated, e.g., = \ (jJc^ +o-^'^)' ^i^d the quantities iF^^g are themselves averages 

ir% = \ 

1 ° 2 

So, the differencing of a TP term will produce two averaged r3r terms. In addition, averaging 
factors before forming products causes the average and difference operators to commute, so that 
the second term in (44) becomes 

1 
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[(.r%-„r^5) +(or«^5--ir«'^5)] 



r«A', 



rP r"'' 



0^ AO' 



5,e 



which exactly cancels the first term on the right side of equation (45). A similar process with 
another — r3r term will cancel the second term on the right side of that equation. 



5.3.4. Non-cancellation of the Terms 

While we are reasonably confident that, with these measures, r3r terms in the Bianchi iden- 
tities will cancel, it is clear that the terms will not cancel to machine accuracy. These are all 
produced by connection of the double-F terms in the Riemann tensor. However, when forming the 
Bianchi identities, the connection takes place on an already-multiplied and summed TT product. 
One cannot undo the sums and products, form averages, and then re-form the triple-F product. And 
the product of averages does not commute with the average of products. The Bianchi identities, 
therefore, will be left with terms proportional to the truncation error and {'dgf' . Unless some clever 
averaging scheme can be found to allow the triple-F terms to also cancel to machine accuracy, any 
CT scheme developed along these lines will be subject to truncation error. The hope, then, is that 
cancellation of second and first order derivatives of the connection (third and second order deriva- 
tives of the metric coefficients) will be sufficient to improve the stability of the discrete evolution 
scheme. 
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Therefore, while the staggered grid CT method clearly works in Riemann normal coordinates, 
detailed numerical experiments will be needed to see if it maintains its desirable properties when 
applied in a general coordinate basis. 



5.4. Numerical Implementation 

5.4.1. General Iterative Approach 

The goal of this CT scheme is to produce an interlocking, staggered four-dimensional grid 
of g(xp tensor values by successively adding new half and whole temporal hypersurfaces. Because 
interpolation is in time, as well as space, the scheme necessarily will be an implicit one and, there- 
fore, iterative. Our approach then will be to first produce an initial guess for g^p on the next level 
of hypersurfaces using an existing explicit scheme. This solution will not satisfy the interlocking 
staggered grid equations, so it will be iterated using the latter until it does. The exact iterative 
scheme is still under development, but one possible implementation is a multigrid technique in 
which the variables are the spatial gij and the equations are the R'^ = field equations. 

This method appears similar to recently- suggested approaches in which the constraints are 
re-solved at each time step. There are, however, two key differences. First, the constraints are not 
solved explicitly. Instead, the evolution equations are solved in such a way that they are implicitly 
enforced. This ensures a scheme in which the evolution and constraints are fully compatible and 

not tracking different numerical solutions to the differential equations. Secondly, the equations 
being iterated are implicit hyperbolic, not elliptic. In addition to information on the new hyper- 
surface, at each iteration these equations utilize information from the previous hypersurface — 
information that constraint equations do not have. With this added information, if implemented 
properly, the iterative scheme should exhibit faster convergence than a regular constraint solver 
would have. 

Figure 5 shows implementation of a simple 1+1 space-time problem. Comparison of this 
figure with Figure 4 will give insight into implementation of the full 4-D interlocking grid, "goo 
and "gii are quantities known from the previous time step, the first being a gauge condition and 
the second being the solution, "''"^gio and "''"^goo are gauge conditions on the new hypersurfaces, 
and "'^^gn is the new solution that will be determined by the iterative numerical scheme. 
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Fig. 5. — Illustration of a 1+1 field evolution, including the initial data problem at and a one- 
point boundary data problem at o^- Bold quantities on the grid are solutions to the initial data and 
evolutionary problems. The computation begins by using gauge conditions to specify goo at ^^t, 
^t, and and goi at and h. Then R^^i, R'^o, and are solved for gn at ^^t, ^t, and ^t. At the 
qx boundary gi 1 may be specified, but at the -ix boundary gi 1 must be consistent with R^o = 0- 
Similarly, the gi 1 at at+ix are obtained by solving R^o = at j^_^_iX. The evolution from to 

proceeds by specifying goi at and goo at gn is computed by solving R^ =0. The 



boundary conditions R^o = also must be applied at ("^^t, 
-ix and yv+i.x: on the new hypersurface. 



x) and ("+2?, ^^i-^) to obtain gn at 
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5.4.2. Using Gauge Conditions and Time Stepping 

Specification of tlie coordinate gauge is similar to doing so in classic 3+1 schemes: goo is 
freely specified at cube corners on whole hypersurfaces, and go, is freely specified at cube edge 
centers on half hypersurfaces. go/ is related to the shift vector (3' = g*^', and both are properly 
located at half time steps vector points, goo = — + g;;P'P-' is related to the lapse a, and properly 
located in this scheme at scalar points. The g/y are the six unique gravitational potential fields that 
we intend to solve. 

A typical time step then begins with g^p components known on hypersurfaces ^~^t,^~^t, and 
One then specifies the go; vector on hypersurface """"^f and the goo scalar at """"^r. The explicit 
predictors for gij at then complete the new g^p field, and the iteration on the gij can then 
begin. If the goa are dependent on the gij field (which generally will be the case here), then the 
gauge conditions need to be updated at each iteration for a consistent solution. 



5.4.3. Constructing the Inverse Metric 

The inverse of gpfp is often needed to raise the connection coefficients for use in the evolution 
equations. In principle these g"P should also be located at second-ranked tensor grid points. In 
practice, however, we never need the actual staggered g"P fields. Instead, we need their interpolated 
averages g"P at third-ranked tensor points. Furthermore, in order for the raising and lowering of 
connection coefficient indices to commute, the staggered field of g"P values must be constructed 
in such a way that its average, and that of g(xp, are orthogonal at those third-ranked tensor points 

t'U = S«P (47) 

So construction of the index-raising tensor g"P is a straightforward matter of averaging g^fp to 
places where the r^py are computed, and then inverting that average locally at those grid points. 
The actual staggered fields of g"P values, whose averages should give us these g"^ at third-ranked 
tensor points, never need to be determined. This procedure is fast, gives us raising and lowering op- 
erators that commute, and follows the afore-mentioned rule of averaging first and then multiplying 
and summing second. 



5.4.4. The Initial Data Problem 

The initial data problem for a staggered grid will be a little more complicated than that for a 
non-staggered scheme. In the latter case, there are 12 unknowns {gij and giy,o) and four equations 
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(R"p = 0) on the initial hypersurface at t = ^t, for a net total of eight degrees of freedom. In the 
staggered case, R^o is computed att= ^t, but the R'^^ are computed att = ~^t. 

Solution of the momentum constraints, given gij at t = ^t, is fairly straightforward in the 
staggered case. They are not functions of second-order time derivatives of the metric (e.g., gij,oo), 
so the placement of those constraints at t = allows them to directly relate at to °g;y. 
All Fs can be computed in a staggered manner. The momentum constraint solution at~^t, then, 
would have six unknowns (^gij) and three constraints, leaving three degrees of freedom, just like 
the non-staggered grid case. 

The Hamiltonian constraint, however, presents a problem. While R^o also does not involve 
any second-order time derivatives, nevertheless it does involve ^rs? order time derivatives gjy,o- 
Those still are defined on half-hypersurfaces and, therefore, need to be interpolated to ^t. That 
is, we need the g, j,o at both "3? and 3f anyway, even if we are not going to compute g/y,oo- At 
first glance there does not appear to be a method of providing an accurate g(y,o(0 fi^^d to properly 
compute this interpolation. Of course, one simply could extrapolate g/y.o at ^^t forward to '^t 
(i.e., assume g/j.oo = 0) and then solve R'^q = there. The Hamiltonian constraint solution then 
would have the six g,y unknowns and one constraint — five degrees of freedom — just like the 
non-staggered case. 

But a serious problem still remains. There is no means of enforcing the Hamiltonian constraint 
att = ^t. While the staggered grid is, in principle, capable of doing that, it can do so only through 
the evolution equations R'y = 0. But at this stage we have not yet begun to enforce the evolution. 
One possible method of solving this is to do nothing. Just accept the fact that, att = ^t, bPq = 
is good only to truncation accuracy. Choosing a very small At (a "thin sandwich") would keep 
this error small. A second approach would be to also explicitly enforce the constraint ont = ^t, 
which would reduce the number of degrees of freedom of the Hamiltonian problem from five to 
four. While producing a more constrained problem, this solution would result in two successive 
hypersurfaces on which the Hamiltonian constraint is satisfied. Yet a third alternative would be to 
locate R'^ j on the 2t hypersurface and extrapolate quantities like gij o forward to ^t. The problem 
with this approach is that, while the Hamiltonian constraint is satisfied for the extrapolated gij,o 
field, when the evolution is begun, the Hamiltonian constraint that is implicit in the evolution 
equations will use fields that are interpolated in time between and ^t. There will be, therefore, 
an implicit constraint violation injected into the evolution at the outset. 

The elegant, and proper, method of solving this problem is to solve all ten of the Einstein 
field equations simultaneously on the initial hypersurfaces. There will be 18 unknowns (gij at ~^t, 
^t, and ^t) and 10 equations, leaving 8 degrees of freedom, just like the non-staggered case. The 
resulting fields will be properly staggered, and the Fs will be properly staggered and interpolated. 
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The Hamiltonian constraint will be satisfied at using the correctly interpolated fields, and it will 
be satisfied at (and even ~^t) as well, because the evolution equations (and therefore the Bianchi 
identities) are fully enforced. Similarly, the momentum constraints will be satisfied at and 
for the same reason, regardless of whether they are actually applied at ^2/^ or at 2 l The reader will, 
of course, recognize that this is more than solving an initial data problem; in actuality the proposed 
scheme solves the initial data problem plus the first file evolutionary time step simultaneously. This 
is done to ensure that the initial data on the first three hypersurfaces are solutions of the discrete, 
staggered evolutionary field equations. No constraint violation will be introduced implicitly other 
than what is naturally present in the evolutionary method already. 



5.4.5. Boundary Conditions 

One-point Boundary Conditions 

When g and n • Vg are specified on the same boundary, where n is the boundary normal, the 
boundary data problem is similar to the initial data problem. In all such cases, the boundary 
constraints are given by 

n^R^p = (48) 

and for rectilinear grids with the boundary normal being a coordinate unit l-form n = w^'^, this 
yields 

R(')p = (49) 

where the symbol (z) is a label indicating the boundary direction in question, not strictly a co- 
ordinate index. The equation R*^'), (no sum) is located at 3-cube corners and plays the role of 
boundary constraint in much the same manner as the Hamiltonian constraint does at Similarly, 
the equations R^'^^ = (a 7^ i) play the same role as the momentum constraints did earlier. By 
analogy, then, the boundary problem is as follows. There are 12 unknowns {gab [a^ i, i] at 
x^') = qx^'^ and atx^'^ = -ix^'\ i.e., on the boundary and one ghost node beyond the boundary), 
and there are 4 equations (49). This leaves 8 degrees of freedom again, which must be specified 
with additional boundary conditions. The diagonal constraint R^'^, = is applied at qx^'^ at 3-cube 
comers. 

Two of the off-diagonal constraints (R^'^; =0 [j ^ i]) are applied on whole hypersurfaces at 

_ix^'\ and the final constraint R^^'^o = is applied also at _ix(') but on time half -hypersurfaces. 

2 2 
The reader will note that this latter set of equations is related to the set of momentum constraints 

r(0);- = 0, but the momentum constraints are defined at ix^'\ 3X^'\ 5X^'\ i^riix^''' and propagated 
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forward by the evolution equations. The constraint at must be expUcitly enforced in most 

cases, along with the three others and the eight freely-specified gap at -ix^'^ and o^^^'^- 

Two-point Boundary Conditions 

Any number of additional boundary data problems are possible. R^^'^a = could be specified at the 
upper (i) boundary, while R('\ = could be be specified at the lower (?) boundary, for example. 
Also, some of the eight free g^p could be specified on opposing boundaries as well. 

For periodic boundary conditions, the six unknown g^^^ at are set to those near nx^'\ 

and those at n+i^^'^ are set to those near ox^'\ No constraints are applied explicitly, only implicitly 
through the wrapping conditions and the evolutionary solution of the field equations. 

5.4.6. Implementation Summary 

It is useful to summarize how a complete problem will proceed. The computation begins by 
using gauge conditions to specify goo at ~^t, ^t, and ^t and go; at ^^t and ^t. The full initial 
data plus time-step problem, including the field equations, is then solved for g, j on ^^t, ^t, and 
^t, applying eight freely-specified gij (or eight functions thereof) in the process. Appropriate 
boundary constraints also need to be applied in order to obtain a consistent solution. 

The field equations are generated as follows. An initial guess for the gij is obtained by some 
means, perhaps using conformal or other existing initial value methods or, for the evolution, an 
explicit forward integration scheme. The full g„p field values then are differenced onto third-rank 
tensor grid points, and the Tf^^y are formed. The g„p values are also averaged to those same third- 
rank tensor points, and an inverse of that average g^^ is used to raise the connection coefficients 
(equations 16 and 25). The Fs, in turn, then are differenced and averaged to second-rank tensor 
points for the Riemann tensor calculation. (Note that there will be no need for values at other 
fourth-rank tensor points [hypercube body centers], as Riemann will be immediately contracted 
into Ricci.) Equation (24) should suffice for the Riemann calculation (i.e., the explicit version 
[equation 41] should not be needed), because our raising and lowering operators commute, ren- 
dering the results of equations (16) and (42) the same to machine accuracy. Contraction to Ricci 
is trivial, as it sums Riemann components already computed at second-rank tensor points, and the 
R'j are then tested to see if they are zero. If not, the local values of gij are modified in an appro- 
priate manner, and the computation of the Ricci components is repeated until convergence. It is 
important to remember that at grid boundaries, the boundary constraints (which are evolutionary 
equations in their own right) must be applied and iterated upon as well, and the freely-specified 
boundary conditions must be applied. 
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When the iterative solution is acceptable, the computation continues to the next hypersurface, 
first specifying the gauge conditions and then solving the field equations. 



6. Discussion 

The largest uncertainty in the proposed method is the effect of the truncation errors intro- 
duced by non-cancellation of the triple-F terms. Analytically we can be sure of exact constraint 
propagation only if the Bianchi identities are satisfied exactly, and that is not the case in a general 
coordinate basis. On the other hand, the fact that the method does work in a local Riemann nor- 
mal system, and has some attractive properties of a finite differential calculus, are encouraging. 
Furthermore, most of the time-dependent derivatives of g do cancel in the general case, and this 
should enhance stability. Analytical investigation of the stability of this method is difficult, so we 
have chosen to do so numerically. 

Two numerical implementations of this scheme are being developed currently. The first as- 
sumes symmetries in two spatial dimensions, rendering the scheme explicit and avoiding the need 
for iteration. The second is a full implementation using staggered grids in four dimensions. Results 
of these studies will help in determining the stability of the method. 

At best, the method is expected to be conditionally stable. The Evans-Hawley and Yee meth- 
ods have this property, and are subject to a Courant-like condition on the time step (De Raedt et al. 2002). 
With c = 1, this implies A? < Ax, which a typical condition applied in most numerical relativity 
implementations. 

If it can be shown that such techniques provide stable constraint transport for finite difference 
implementations of numerical relativity, then other implementations might benefit from similar 
schemes that use four-dimensional grids and that have temporal and spatial derivatives that com- 
mute. In the finite element case, this would begin by extending the elements into the time direction, 
rather than simply time- stepping a three-dimensional finite element grid. However, some addi- 
tional features would have to be introduced, analogous to grid staggering, to ensure commutation 
of spatial and time derivatives across element boundaries. 

The pseudo-spectral may be intractable at the present time. As spatial derivatives are com- 
puted using the full spatial extent of the grid, in order to create a method in which these commuted 
with time derivatives, one could imagine using many, if not all, previous and future time hyper- 
surfaces to compute the latter. This would be a truly four-dimensional grid method and involve 
solving the entire spacetime structure in one giant iterative procedure. Present-day computers still 
struggle with three-dimensional explicit schemes, so a four-dimensional implicit one is clearly 
beyond current technology. However, in the not-too-distant future such codes might begin to be 
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feasible. 

The author has benefited greatly from discussions with F. Estabrook, L. Lindblom, M. Miller, 
D. Murphy, and M. Scheel. He is also grateful for a JPL Institutional Research and Development 
grant, and for the continued hospitality of the TAPIR group at Caltech. This research was per- 
formed at the Jet Propulsion Laboratory, California Institute of Technology, under contract to the 
National Aeronautics and Space Administration. 

REFERENCES 

Clarke, D. A. 1996, ApJ, 457, 291-320. 

De Raedt, H., Kole, J. S., Michielsen, K. F. L., & Figge, M. T. 2002, arXiv:physics/02 10035. 
Evans, C. R. & Hawley, J. F 1988, ApJ, 332, 659-677. 
Stone, J. M. & Norman, M. L. 1992a, ApJS, 80, 753-790. 
Stone, J. M. & Norman, M. L. 1992b, ApJS, 80, 791-818. 

Scheel, M. A., Kidder, L. E., Lindblom, L., Pfeiffer, H. P, & Teukolsky, S. A. 2002, Phys. Rev. D, 
66, 124005. 

Yee, K. S. 1966, IEEE Transactions on Antennas and Propagation, 14, 302-307. 



This preprint was prepared with the AAS MTeX macros v5.0. 



